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Abstract 


In a recent paper in Astrophysics and Space Science Vol. 364 no. 11 (2019), S. Ershkov & D. 
Leschenko presented a new solving procedure for Euler-Poisson equations for solving momentum 
equations of the CR3BP near libration points for uniformly rotating planets having inclined orbits in 
the solar system with respect to the orbit of the Earth. The system of equations of the CR3BP has been 
explored with regard to the existence of an analytic way of presentation of the approximated solution 
in the vicinity of libration points. A new and elegant ansatz has been suggested in their publication 
whereby, in solving, the momentum equation is reduced to a system two coupled Riccati ODEs. In 
this paper, we presented a numerical solution of such coupled Riccati ODEs using Mathematica 
software package. 
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1. Introduction. 


The equations of motion of the R3BP describe the dynamics of infinitesimal mass 
m under the action of gravitational forces affected by two celestial bodies of giant 
masses (in this problem M sun and m planet, m planet < M Sun), which are rotating around 
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their common centre of masses on Kepler’s trajectories. The aforesaid infinitesimal 
mass m is supposed to be moving (as first approximation) inside the restricted region 
of space around the planet of mass m planet or inside the so-called Hill sphere (where 


dp iS Semimajor axis of the planet). [1] 


It is worth to note that a large amount of previous and recent works concerning 
analytical development with respect to these equations exist. We can mention a few 
of them here. 

Motivated by trajectory design applications, analytical work from the 1970s 
includes Farquhar and Kamel’s third-order expansion of orbits in the vicinity of the 
Earth— Moon L2 libration point using a Poincar’e-Lindstedt method.5 During the 
same time period, Richardson and Cary developed a third-order approximation of 
quasi-periodic motion near the Sun—Earth L1 and L2 libration points via the method 
of multiple time scales. More recent semi-analytical work leverages computer algebra 
tools. G’omez et al. use a Poincar’eLindstedt method to generate high-order quasi- 
halo orbit expansions. Jorba and Masdemont investigate the dynamics around the 
libration points using center manifold reduction. These local methods offer a thorough 
view of the dynamics in the libration point vicinity, but they are limited by their 
regions of convergence. Specialized algebraic manipulators are often required as 
well, which can be an implementation obstacle. Still, the solutions have been 
successfully exploited for generating heteroclinic connections.[5] 

In the 1980s Howell and Pernicka presented a numerical shooting approach for 
correcting approximate quasi-periodic orbits, though lacking control of orbit 
parameters. Gomez and Mondelo later developed a scheme for computing two- 
dimensional quasi-periodic tori by using an invariant circle parameterized by Fourier 
coefficients of a stroboscopic map. Kolemen et al. use a similar approach except with 
multiple Poincare maps and directly parameterizing the invariant circle using states. | 
The current approach uses similar concepts to these methods but with a more general 
formulation. We use a stroboscopic map, which requires less knowledge of the 
solution structure than a sequence of Poincare maps, and consider tori of arbitrary 
dimension. We use a grid of states on the invariant torus allowing us to avoid a 
convolution that is necessary when using a Fourier coefficient parameterization. In 
addition, we incorporate a different set of solution and continuation constraints that 


are broadly applicable. This allows the method to be used unchanged between 
2 


families of solutions. An alternative approach presented by Schilder et a/. computes 
a torus of a flow directly, and some of the constraints used have been adapted from 
this approach. The flow approach has been applied to the CR3BP; however, this 
requires dealing with a torus of dimension one larger than the torus of an associated 


map.[5] 


We should especially emphasize the theory of orbits, which was developed in 
profound work [4] for the case of the Circular Restricted Problem of Three Bodies 
(CR3BP) (primaries M sun and m planet are rotating around their common centre of 
masses on circular orbits). According to [4], equations of motion of the infinitesimal 
mass m should be presented in the co-rotating frame of a Cartesian coordinate system 
7 ={x, y, z} in case of CR3BP (at given initial conditions). 

The suggested ansatz could be useful from a practical point of view in celestial 
mechanics (for example, to obtain a regular data of astrometric observations). Indeed, 
sometimes it is convenient for an observer to consider such the celestial dynamical 
system at a fixed angle to the plane of mutual circular rotation of the chosen primaries 
e.g. in the dynamics of planetary rotation in the Solar system.[1] 

The most common example is Pluto’s angle of orbit’s inclination with respect to 
Earth’s orbit (which is 17°9’); namely, Pluto's orbit is inclined, or tilted, circa 17.1 
degrees from the ecliptic of Earth’s orbit. Except Mercury's inclination of 7 degrees, 


all the other orbits of planets are closer to the orbit of Earth.[1] 


2. Coupled Riccati ODEs. 


According to the ansatz for solving Poisson equations, the system of Equations for 
three problem has the analytical way to present the general solution in regard to the 


time f:[1] 
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here oO = O(x, y, Z) is some arbitrary (real) function, given by the initial conditions; 
the real-valued coefficients a(x,y,z, ft), b(x,y,z, t) in (1) are solutions of the mutual 
system of two Riccati ordinary differential equations, ODEs in regard to the time ¢ 


(here & = 2): 
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System (2) describes the evolution of function a in dependence on the function b 
in regard to the time ¢ (and vice versa); we should especially note that the Riccati 


ODE has no analytical solution in general case. 


3. Numerical solution of equation (2) 


As an alternative approach to aforementioned solution as described in [1], we 
would like to prove that numerical solutions exist for equation (2). 


We will obtain the numerical solution using Mathematica software package. 


Case 1: 


B=k-w,=0.25 
g=k-w,=0.15 
ybar =k-w, =0.25 


param={B->0.15,g->0.15,ybar->0.25,y->1.0,5->0.03,Pstar->0.025, iM->0}; 

eqn1=P'[t]==(B/2)*P[t]*2-(g*c[t])*P[t]-(B/2)*(c[t]*2-1)+ybar*c[t]/.param; 
eqn2=c'[t]==(-(g/2)*c[t]*2+(B*P[t])*c[t]+(g/2)*(P[t]*2-1)- 

ybar*P[t])/.param; 

sol=Flatten[NDSolve[{eqn1,eqn2, P[O]==-0.15,c[0]==0.7},{P, c},{t,-20,10}]]; 

plt1=ParametricPlot[{P[t],c[t]}/.sol,{t, -7,7}, Frame->True, FrameLabel- 


>{Style["m",18],Style["c", 18]},AspectRatio->1] 


sol=Flatten[NDSolve[{eqni1,eqn2, P[O0]==-0.1,c[0]==0.7},{P, c},{t,-20,10}]]; 
plt2=ParametricPlot[{P[t],c[t]}/.sol,{t, -7,7}, Frame->True, FrameLabel- 
>{Style["r",18],Style["c", 18]},AspectRatio->1]; 
sol=Flatten[NDSolve[{eqn1,eqn2, P[0]==-.07,c[0]==0.7},{P, c},{t,-30,15}]] 


plt3=ParametricPlot[{P[t],c[t]}/.sol,{t, -10,10}, Frame->True, FrameLabel- 
>{Style["m",18],Style["c", 18]},AspectRatio->1]; 
sol=Flatten[NDSolve[{eqn1,eqn2, P[0]==-.06,c[0]==0.7},{P, c},{t,-30,15}]] 


plt4=ParametricPlot[{P[t],c[t]}/.sol,{t, -10,10}, Frame->True, FrameLabel- 
>{Style["r",18],Style["c", 18]}, AspectRatio->1]; 
sol=Flatten[NDSolve[{eqn1,eqn2, P[O0]==-.05,c[0]==0.7},{P, c},{t,-30,15}]] 
plt5=ParametricPlot[{P[t],c[t]}/.sol,{t, -10,10}, Frame->True, FrameLabel- 
>{Style["m",18],Style["c", 18]},AspectRatio->1]; 


Show[{plt1, plt2, plt3, plt4, plt5}, AspectRatio > 1, Axes > False] 
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Case 2: 
B=k-w, =0.75 
g=k-w,=0.35 


ybar =k-w, =0.25 


param={B->0.75,g->0.35,ybar->0.25,y->1.0,5->0.03,Pstar->0.025, iM->0}; 

eqn1=P'[t]==(8/2)*P[t]*2-(g*c[t])*P[t]-(B/2)*(c[t]*2-1)+ybar*c[t]/.param; 
eqn2=c'[t]==(-(g/2)*c[t]*2+(B*P[t])*c[t]+(g/2)*(P[t]*2-1)- 

ybar*P[t])/.param; 

sol=Flatten[NDSolve[{eqn1,eqn2, P[O]==-0.15,c[0]==0.7},{P, c},{t,-20,10}]]; 

plt1=ParametricPlot[{P[t],c[t]}/.sol,{t, -7,7}, Frame->True, FrameLabel- 


>{Style["m",18],Style["c", 18]},AspectRatio->1] 


sol=Flatten[NDSolve[{eqni1,eqn2, P[0]==-0.1,c[0]==0.7},{P, c},{t,-20,10}]]; 
plt2=ParametricPlot[{P[t],c[t]}/.sol,{t, -7,7}, Frame->True, FrameLabel- 
>{Style["r",18],Style["c", 18]}, AspectRatio->1]; 
sol=Flatten[NDSolve[{eqni1,eqn2, P[0]==-.07,c[0]==0.7},{P, c},{t,-30,15}]] 


plt3=ParametricPlot[{P[t],c[t]}/.sol,{t, -10,10}, Frame->True, FrameLabel- 
>{Style["m",18],Style["c", 18]}, AspectRatio->1]; 
sol=Flatten[NDSolve[{eqni1,eqn2, P[0]==-.06,c[0]==0.7},{P, c},{t,-30,15}]] 
plt4=ParametricPlot[{P[t],c[t]}/.sol,{t, -10,10}, Frame->True, FrameLabel- 
>{Style["r",18],Style["c", 18]},AspectRatio->1]; 
sol=Flatten[NDSolve[{eqn1,eqn2, P[O0]==-.05,c[0]==0.7},{P, c},{t,-30,15}]] 
plt5=ParametricPlot[{P[t],c[t]}/.sol,{t, -10,10}, Frame->True, FrameLabel- 
>{Style["r",18],Style["c", 18]},AspectRatio->1]; 


Show[{plt1, plt2, plt3, plt4, plt5}, AspectRatio > 1, Axes > False] 
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Case 3: 
B=k-w,=0.15 
g=k-w,=0.15 


ybar =k-w, =0.01 


param={B->0.15,g->0.15,ybar->0.01,y->1.0,5->0.03,Pstar->0.025, iM->0}; 

eqn1=P'[t]==(8/2)*P[t]42-(g*c[t])*P[t]-(B/2)*(c[t]*2-1)+ybar*c[t]/.param; 
eqn2=c'[t]==(-(g/2)*c[t]*2+(B*P[t])*c[t]+(g/2)*(P[t]}*2-1)- 

ybar*P[t])/.param; 

sol=Flatten[NDSolve[{eqn1,eqn2, P[O]==-0.15,c[0]==0.7},{P, c},{t,-20,10}]]; 

plt1=ParametricPlot[{P[t],c[t]}/.sol,{t, -7,7}, Frame->True, FrameLabel- 


>{Style["m",18],Style["c", 18]},AspectRatio->1] 
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sol=Flatten[NDSolve[{eqni1,eqn2, P[0]==-0.1,c[0]==0.7},{P, 
c},{t,-20,10}]]; 
plt2=ParametricPlot[{P[t],c[t]}/.sol,{t, -7,7}, Frame->True, FrameLabel- 
>{Style["r",18],Style["c", 18]},AspectRatio->1]; 
sol=Flatten[NDSolve[{eqni1,eqn2, P[O]==-.07,c[0]==0.7},{P, c},{t,-30,15}]] 


plt3=ParametricPlot[{P[t],c[t]}/.sol,{t, -10,10}, Frame->True, FrameLabel- 
>{Style["m",18],Style["c", 18]},AspectRatio->1]; 
sol=Flatten[NDSolve[{eqni1,eqn2, P[0]==-.06,c[0]==0.7},{P, c},{t,-30,15}]] 
plt4=ParametricPlot[{P[t],c[t]}/.sol,{t, -10,10}, Frame->True, FrameLabel- 
>{Style["r",18],Style["c", 18]}, AspectRatio->1]; 
sol=Flatten[NDSolve[{eqn1,eqn2, P[O0]==-.05,c[0]==0.7},{P, c},{t,-30,15}]] 
plt5=ParametricPlot[{P[t],c[t]}/.sol,{t, -10,10}, Frame->True, FrameLabel- 
>{Style["m",18],Style["c", 18]},AspectRatio->1]; 


Show[{plt1, plt2, plt3, plt4, plt5}, AspectRatio > 1, Axes — False] 


The above shematic diagrams show numerical results of equation (2) in three different 


sets of parameters. 


4. Discussion and conclusion. 


As we can see in [1], the system of equations for the Circular Restricted Problem 
of Three Bodies (CR3BP) (which governs the dynamics of infinitesimal mass m under 
the action of gravitational forces affected by two mutually rotating celestial bodies of 
giant masses) has been proved to be very complicated to solve analytically in 3D case. 

Indeed, numerical solutions of the mutual system (2) of two Riccati ordinary 
differential equations in regard to the time ¢, can be found using Mathematica. 
We also note that due to the special character of the solutions of Riccati-type ODEs, 
there is the possibility for sudden jumping in the magnitude of the solution at some 
time fo (see e.g. the solutions of Riccati-type for Abel ODE of 1-st order). In the 
physical sense, such jumping of Riccati-type solutions could be associated with the 


effect of a sudden acceleration/deceleration of the celestial body’s velocity at a 
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definite moment of parametric time fo. 
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